Ramalho et al 2021
Usefull librarys.
library(rmarkdown) # You need this library to run this template.
library(epuRate) # Install with devtools: install_github("holtzy/epuRate", force=TRUE)
library(DT) # pour faire de joies affichage de dataframe
library(readtext)
#library(tcltk) # choose directory
library(readxl)
library(openxlsx)
library(ggplot2)
library(tidyverse)
# detach("package:MASS")
library(afex) # nouveau package, surveiller les MAJ
library(lme4)
library(RVAideMemoire)
library(LMERConvenienceFunctions)
library(phia)
library(multcomp)
library(lmerTest)
library(emmeans)
library(ggpubr)
library(ggeffects)
library(splines)
library(DHARMa)
library(glmmTMB)
library(sjPlot)
library(sjmisc)
library(rstanarm)
library(brms) # for models
library(bayesplot)
library(bayestestR)Ramalho B. L., Moly J., Raffin E., Bouet R., Harquel S., Farnè A., Reilly K.T. (2021) Electrocutaneous stimulation to the face inhibits motor evoked potentials in the hand: Face-hand sensorimotor interactions revealed by afferent inhibition. EJN
https://pubmed.ncbi.nlm.nih.gov/34796553/
In this study we used the afferent inhibition protocol to assess sensorimotor interactions between the hand and body parts represented close to the hand area in the sensorimotor cortex: the face and the upper limb. We found that the amplitude of MEPs in the right FDI was reduced when the TMS pulse was preceded by an electrocutaneous stimulus on the face (upper lip or cheek) or on the upper limb (arm or forearm), but this inhibition was more robust when associated with face stimulation. These results provide the first evidence for afferent inhibition between the skin on the face and a hand muscle, but also between the skin on the upper limb and a hand muscle.
Key words:
Time series; GLM; Gamma; Splines; Power; excel sheets;
Datas
Extract
Here we define a function to extract datas from exel files:
- Read excel file
- Explore all sheet (coresponding to subjects)
- TP stimuation become O ms timing.
Like this I can compare each time point (15,15,35….) with 0 and we keep
TP variance.
- Each mesure was centered (\(Mesure.TP =
Mesure / TP\))
Exitability difference with Test Pulse baseline.
We add one factor : Localization
LIP & CHEEK : loc1
FOREAM & ARM : loc2
read.all_subject <- function(prefix){
filename <- paste("/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Datas_raw/", prefix, "_all_subjects.xlsx", sep = "")
sheets.names <- excel_sheets(filename)
datas.sheets <- lapply(excel_sheets(filename), read_excel, path = filename)
datas.raw <- data.frame()
# first sheet is MEAN
for (xi_suj in 1:length(datas.sheets)){
datas.raw <- bind_rows(datas.raw, mutate(datas.sheets[[xi_suj]], Subject = sheets.names[xi_suj]) )
}
rm(xi_suj)
rm(filename, sheets.names, datas.sheets)
# Extract subject TP mean
datas.raw %>%
dplyr::select(Subject, TP) %>%
group_by(Subject) %>%
summarise(TP_mean = mean(TP),
TP_sd = sd(TP),
TP_med = median(TP)) -> datas.TP
datas.raw %>%
mutate("0" = TP) %>%
dplyr::select(-TP) %>%
gather(key = "Time", value = "Mesure", -Subject) %>%
mutate(Time = as.factor(Time)) %>%
drop_na() %>%
left_join(datas.TP) %>%
group_by(Subject, Time, Mesure) %>%
data.frame() %>%
mutate(Mesure.TP = Mesure/TP_mean,
Mesure.TP.med = Mesure/TP_med) %>%
mutate(Body = as.factor(prefix)) %>%
mutate(Subject = as.factor(Subject)) -> datas.raw
return(datas.raw)
}
datas.forearm <- read.all_subject("FOREARM")
datas.arm <- read.all_subject("ARM")
datas.cheek <- read.all_subject("CHEEK")
datas.lip <- read.all_subject("LIP")
bind_rows(datas.forearm, datas.arm, datas.cheek, datas.lip) %>%
mutate(Subject = as.factor(Subject),
Body = as.factor(Body)) %>%
dplyr::select(Body, Subject, Time, Mesure, Mesure.TP, Mesure.TP.med) -> datas.raw
# Add localization factor
datas.raw %>%
mutate(Localization = case_when(Body == "LIP" ~ "Loc1",
Body == "CHEEK" ~ "Loc1",
Body == "FOREARM" ~ "Loc2",
Body == "ARM" ~ "Loc2")) %>%
mutate(Localization = as.factor(Localization)) -> datas.raw
datas.raw %>%
dplyr::select(-Subject) %>%
gtsummary::tbl_summary(by = "Body") %>%
gtsummary::add_p()| Characteristic | FOREARM, N = 2,0141 | ARM, N = 1,7201 | CHEEK, N = 1,5181 | LIP, N = 1,8241 | p-value2 |
|---|---|---|---|---|---|
| Time | >0.9 | ||||
| 0 | 476 (24%) | 410 (24%) | 346 (23%) | 416 (23%) | |
| 15 | 195 (9.7%) | 166 (9.7%) | 147 (9.7%) | 178 (9.8%) | |
| 25 | 190 (9.4%) | 166 (9.7%) | 148 (9.7%) | 179 (9.8%) | |
| 35 | 185 (9.2%) | 163 (9.5%) | 152 (10%) | 180 (9.9%) | |
| 45 | 193 (9.6%) | 159 (9.2%) | 143 (9.4%) | 172 (9.4%) | |
| 55 | 193 (9.6%) | 165 (9.6%) | 144 (9.5%) | 167 (9.2%) | |
| 65 | 190 (9.4%) | 163 (9.5%) | 147 (9.7%) | 172 (9.4%) | |
| 75 | 199 (9.9%) | 165 (9.6%) | 142 (9.4%) | 180 (9.9%) | |
| 85 | 193 (9.6%) | 163 (9.5%) | 149 (9.8%) | 180 (9.9%) | |
| Mesure | 0.78 (0.39, 1.37) | 0.87 (0.46, 1.36) | 0.58 (0.24, 1.17) | 0.58 (0.27, 1.13) | <0.001 |
| Mesure.TP | 0.81 (0.41, 1.23) | 0.88 (0.46, 1.31) | 0.72 (0.26, 1.36) | 0.72 (0.35, 1.35) | <0.001 |
| Mesure.TP.med | 0.90 (0.49, 1.38) | 0.90 (0.48, 1.33) | 0.83 (0.32, 1.55) | 0.94 (0.51, 1.75) | <0.001 |
| Localization | <0.001 | ||||
| Loc1 | 0 (0%) | 0 (0%) | 1,518 (100%) | 1,824 (100%) | |
| Loc2 | 2,014 (100%) | 1,720 (100%) | 0 (0%) | 0 (0%) | |
| 1 n (%); Median (IQR) | |||||
| 2 Pearson's Chi-squared test; Kruskal-Wallis rank sum test | |||||
rm(datas.arm, datas.cheek, datas.forearm, datas.lip, read.all_subject)Display
Distribution
Some displays to explore distributions. Here time is categorical factor.
datas.raw %>%
mutate(Time = Time) %>%
ggplot(aes(x = Time, y = Mesure.TP, color = Time)) +
geom_jitter(color = "black", alpha = 0.1) +
geom_violin(alpha = 0, position = position_dodge(width = 0.9)) +
geom_boxplot(width=.2, alpha = 0, outlier.alpha = 0,
position = position_dodge(width = 0.9)) +
stat_summary(fun.y = "mean",
geom = "point",
size = 2,
position = position_dodge(width=0.9),
color = "black") +
coord_cartesian(ylim =c(-0.1, 5)) +
facet_wrap(~Body)Datas does’t looks Gaussian.
Focus on “ARM” :
datas.raw %>%
mutate(Time = Time) %>%
filter(Body == "ARM") %>%
ggplot(aes(x = Time, y = Mesure.TP, color = Time)) +
geom_violin(alpha = 0, position = position_dodge(width = 0.9)) +
geom_boxplot(width=.2, alpha = 0, outlier.alpha = 0,
position = position_dodge(width = 0.9)) +
stat_summary(fun.y = "mean",
geom = "point",
size = 2,
position = position_dodge(width=0.9),
color = "black") +
coord_cartesian(ylim =c(-0.1, 4)) +
facet_wrap(~Subject)
As we can see. It’s usefull to add a “0” time as a TP periode. Variance
within TP condition help to see a difference between another time.
Variability
Explore inter-subject variability
datas.sum <- Rmisc::summarySE(datas.raw, measurevar="Mesure.TP", groupvars=c("Time","Body"))
pd <- position_dodge(0.3)
ggplot(datas.sum, aes(x = Time, y = Mesure.TP, colour = Body, group = Body)) +
geom_errorbar(aes(ymin = Mesure.TP-ci, ymax = Mesure.TP+ci),
colour = "black", alpha = 0.4,
width = .1, position = pd) +
geom_line(position = pd) +
geom_point(position = pd, size = 3) +
ggtitle("Subject average")datas.raw %>%
group_by(Subject, Body, Time) %>%
summarise(Mesure.TP = mean(Mesure.TP)) %>%
mutate(Time = as.numeric(as.character(Time))) %>%
ggplot(aes(x = Time, y = Mesure.TP, color = Subject)) +
geom_line(width=.2, show.legend = FALSE) +
facet_wrap(~Body, scales = "free") +
theme(legend.position="bottom")+
ggtitle("Subject datas")High Subject’s heterogeneity for Mesure.TP along Time.
Stats
I did 2 analysis:
- like the Karen’s paper, I use Time as a categorial variable
- To introduce a inter-dependance between time periode, also I use Time as a numerical variable
We use a glm method because we have to:
- Consider the Gamma distribution (Datas are continuous and bounded) -
We accound of inter-Subject variability and Time by Body variability
(random effect).
Time as factor
First, Time is consider as factor.
So we explore 2 fixe effects and there interation:
- Body (4 levels)
- Time (9 levels)
We consider 3 random effets:
- Subject intercept
- Body intercept - Time slope by Body
Models
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Fact.RData")
# model.03 <- glmer(Mesure.TP ~ (Body + Time)^2 + (1|Subject) + (Time|Body),
# family = Gamma("identity"),
# data = datas.raw)
# model.03 <- glmer(Mesure.TP ~ (Body + Time)^2 + (1|Subject) + (Time|Body) + (1|Localization),
# family = Gamma("identity"),
# data = datas.raw)
#save(model.01, model.02, model.03,
# file = "/Volumes/crnldata/dycog/Epilepto/Karen/TMS/Scripts/Models_Time_Fact.RData")
Anova(model.03, type = "II")## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Mesure.TP
## Chisq Df Pr(>Chisq)
## Body 10.902 3 0.012270 *
## Time 132.020 8 < 2.2e-16 ***
## Body:Time 57.385 24 0.000148 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We have a significative interaction Body:Time.
Post-hoc
Effect’s displays for raw datas and estimate mesures.
datas.sum <- Rmisc::summarySE(datas.raw, measurevar="Mesure.TP", groupvars=c("Time","Body"))
pd <- position_dodge(0.3)
ggplot(datas.sum, aes(x = Time, y = Mesure.TP, colour = Body, group = Body)) +
geom_errorbar(aes(ymin = Mesure.TP-ci, ymax = Mesure.TP+ci),
colour = "black", alpha = 0.4,
width = .1, position = pd) +
geom_line(position = pd) +
geom_point(position = pd, size = 3) +
ggtitle("Raw datas")plot_model(model.03,
type = "pred",
terms = c("Time", "Body")) +
ggtitle("Estimate datas")Here time-pairwise by Body, Bonferroni correction:
em <- emmeans(model.03, ~ Time | Body)
#CLD(em)
# contrast(em, interaction = c( "pairwise"), adjust = "bonferroni")
des_cont <- list("0_15" = c(-1, 1, 0, 0, 0, 0, 0, 0, 0),
"0_25" = c(-1, 0, 1, 0, 0, 0, 0, 0, 0),
"0_35" = c(-1, 0, 0, 1, 0, 0, 0, 0, 0),
"0_45" = c(-1, 0, 0, 0, 1, 0, 0, 0, 0),
"0_55" = c(-1, 0, 0, 0, 0, 1, 0, 0, 0),
"0_65" = c(-1, 0, 0, 0, 0, 0, 1, 0, 0),
"0_75" = c(-1, 0, 0, 0, 0, 0, 0, 1, 0),
"0_85" = c(-1, 0, 0, 0, 0, 0, 0, 0, 1))
contrast(em, des_cont, adjust = "bonferroni")## Body = FOREARM:
## contrast estimate SE df z.ratio p.value
## 0_15 0.0454 0.0786 Inf 0.577 1.0000
## 0_25 -0.0324 0.0745 Inf -0.435 1.0000
## 0_35 -0.1218 0.0694 Inf -1.754 0.6359
## 0_45 -0.1484 0.0686 Inf -2.164 0.2435
## 0_55 -0.1356 0.0692 Inf -1.959 0.4012
## 0_65 -0.1845 0.0668 Inf -2.761 0.0460
## 0_75 -0.0254 0.0747 Inf -0.340 1.0000
## 0_85 -0.0760 0.0720 Inf -1.055 1.0000
##
## Body = ARM:
## contrast estimate SE df z.ratio p.value
## 0_15 0.1280 0.0861 Inf 1.486 1.0000
## 0_25 -0.0933 0.0738 Inf -1.264 1.0000
## 0_35 -0.1907 0.0706 Inf -2.700 0.0555
## 0_45 -0.3533 0.0621 Inf -5.686 <.0001
## 0_55 -0.4000 0.0636 Inf -6.293 <.0001
## 0_65 -0.2764 0.0715 Inf -3.863 0.0009
## 0_75 -0.1191 0.0786 Inf -1.516 1.0000
## 0_85 -0.0596 0.0824 Inf -0.723 1.0000
##
## Body = CHEEK:
## contrast estimate SE df z.ratio p.value
## 0_15 -0.0981 0.0651 Inf -1.507 1.0000
## 0_25 -0.1923 0.0628 Inf -3.063 0.0175
## 0_35 -0.1439 0.0650 Inf -2.214 0.2145
## 0_45 -0.0639 0.0690 Inf -0.926 1.0000
## 0_55 -0.1945 0.0617 Inf -3.151 0.0130
## 0_65 -0.1228 0.0653 Inf -1.881 0.4803
## 0_75 -0.1221 0.0655 Inf -1.865 0.4972
## 0_85 -0.1843 0.0630 Inf -2.927 0.0274
##
## Body = LIP:
## contrast estimate SE df z.ratio p.value
## 0_15 0.0164 0.0736 Inf 0.222 1.0000
## 0_25 -0.0379 0.0724 Inf -0.523 1.0000
## 0_35 -0.0592 0.0723 Inf -0.819 1.0000
## 0_45 -0.2962 0.0603 Inf -4.913 <.0001
## 0_55 -0.3599 0.0612 Inf -5.885 <.0001
## 0_65 -0.3320 0.0593 Inf -5.596 <.0001
## 0_75 -0.1459 0.0660 Inf -2.210 0.2167
## 0_85 -0.0735 0.0700 Inf -1.051 1.0000
##
## P value adjustment: bonferroni method for 8 tests
First analysis found significant :
- LIP inhibition : 45, 55 and
65 ms
- CHEEK inhibition : 25, 55 and 85 ms
- FOREARM 65 ms
- ARM 45, 55 and 65 ms
Time as numeric
Is intuitive to consider a Time as numeric values,
not factor. There is a dependence between levels.
Also, Time is not a linear effect. So we fit a
non-linear Time effect with spline with 7
breakpoints.
Models
We use exactly the same model as previously.
datas.raw %>%
mutate(Time = as.numeric(as.character(Time))) -> datas.raw.num
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume.RData")
# model.num.04 <- mixed(Mesure.TP ~ (Body + ns(Time,7))^2 + (1|Subject) + (ns(Time,7)|Body),
# family = Gamma("identity"),
# data = datas.raw.num,
# method = "LRT")
# save(model.num.04,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume.RData")
Anova(model.num.04$full_model, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Mesure.TP
## Chisq Df Pr(>Chisq)
## (Intercept) 8.4606 1 0.003629 **
## Body 1.5617 3 0.668099
## ns(Time, 7) 122.1927 6 < 2.2e-16 ***
## Body:ns(Time, 7) 54.1249 18 1.755e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model.num.04)## Mixed Model Anova Table (Type 3 tests, LRT-method)
##
## Model: Mesure.TP ~ (Body + ns(Time, 7))^2 + (1 | Subject) + (ns(Time,
## Model: 7) | Body)
## Data: datas.raw.num
## Df full model: 66
## Df Chisq Chi Df Pr(>Chisq)
## Body 66 0.0052 0
## ns(Time, 7) 60 16.5034 6 0.01129 *
## Body:ns(Time, 7) 48 19.2304 18 0.37776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We always find a significant interaction Body:Time
Post-hoc
Effect’s displays for raw datas and estimate mesures.
datas.sum <- Rmisc::summarySE(datas.raw, measurevar="Mesure.TP", groupvars=c("Time","Body"))
pd <- position_dodge(0.3)
ggplot(datas.sum, aes(x = Time, y = Mesure.TP, colour = Body, group = Body)) +
geom_smooth(width=.2, show.legend = TRUE, se = FALSE) +
ggtitle("Raw datas")plot_model(model.num.04$full_model,
type = "pred",
terms = c("Time", "Body")) +
ggtitle("Estimate datas")There is an effect, ok, but where ? In this case, When ?
To explore this issue, we estimate Mesure.TP for
differents Time periode ((0, 15, 25, 35, 45, 55, 65,
75, 85). In this way, as a factor analysis, we compare each 8 time
periode to 0 (baseline). Like this, we could express when
Mesure.TP is significantly different from baseline.
em <- emmeans(model.num.04, ~ Time | Body,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)),
type = "pred")
des_cont <- list("0_15" = c(1, -1, 0, 0, 0, 0, 0, 0, 0),
"0_25" = c(1, 0, -1, 0, 0, 0, 0, 0, 0),
"0_35" = c(1, 0, 0, -1, 0, 0, 0, 0, 0),
"0_45" = c(1, 0, 0, 0, -1, 0, 0, 0, 0),
"0_55" = c(1, 0, 0, 0, 0, -1, 0, 0, 0),
"0_65" = c(1, 0, 0, 0, 0, 0, -1, 0, 0),
"0_75" = c(1, 0, 0, 0, 0, 0, 0, -1, 0),
"0_85" = c(1, 0, 0, 0, 0, 0, 0, 0, -1))
contrast(em, des_cont, adjust = "bonferroni")## Body = FOREARM:
## contrast estimate SE df z.ratio p.value
## 0_15 0.11936 0.0592 Inf 2.016 0.3507
## 0_25 0.16543 0.0543 Inf 3.049 0.0184
## 0_35 0.15410 0.0569 Inf 2.707 0.0543
## 0_45 0.10241 0.0599 Inf 1.710 0.6982
## 0_55 0.14039 0.0547 Inf 2.568 0.0819
## 0_65 0.16410 0.0571 Inf 2.872 0.0326
## 0_75 0.10402 0.0663 Inf 1.569 0.9340
## 0_85 0.18853 0.0628 Inf 3.004 0.0213
##
## Body = ARM:
## contrast estimate SE df z.ratio p.value
## 0_15 -0.05211 0.0721 Inf -0.723 1.0000
## 0_25 0.04108 0.0603 Inf 0.681 1.0000
## 0_35 0.12770 0.0621 Inf 2.058 0.3171
## 0_45 0.12296 0.0644 Inf 1.910 0.4493
## 0_55 0.16844 0.0572 Inf 2.943 0.0260
## 0_65 0.16148 0.0630 Inf 2.563 0.0830
## 0_75 0.03861 0.0723 Inf 0.534 1.0000
## 0_85 0.07519 0.0722 Inf 1.041 1.0000
##
## Body = CHEEK:
## contrast estimate SE df z.ratio p.value
## 0_15 -0.08051 0.0739 Inf -1.089 1.0000
## 0_25 0.03605 0.0635 Inf 0.568 1.0000
## 0_35 0.21975 0.0626 Inf 3.508 0.0036
## 0_45 0.35782 0.0588 Inf 6.081 <.0001
## 0_55 0.38872 0.0571 Inf 6.806 <.0001
## 0_65 0.29071 0.0648 Inf 4.483 0.0001
## 0_75 0.11239 0.0778 Inf 1.444 1.0000
## 0_85 0.05736 0.0829 Inf 0.692 1.0000
##
## Body = LIP:
## contrast estimate SE df z.ratio p.value
## 0_15 0.00231 0.0667 Inf 0.035 1.0000
## 0_25 0.00442 0.0595 Inf 0.074 1.0000
## 0_35 0.08749 0.0624 Inf 1.402 1.0000
## 0_45 0.27856 0.0567 Inf 4.912 <.0001
## 0_55 0.37798 0.0523 Inf 7.223 <.0001
## 0_65 0.32072 0.0560 Inf 5.727 <.0001
## 0_75 0.15236 0.0643 Inf 2.368 0.1430
## 0_85 0.07082 0.0703 Inf 1.007 1.0000
##
## P value adjustment: bonferroni method for 8 tests
For instance, here we see a significative decrease for LIP at 45, 55 and 65 (p.value was control for multiple comparaison by bonferroni method).
Supplementary analysis
Following the advice of the expert reviewers from this first round of reviews we made substantial changes. All these analysis was recommended by reviewers.
Experiment Factor
We keep model with Time as numeric.
We add factor to code the experiment factor.
LIP & CHEEK : Experiment 1
FOREAM & ARM : Experiment 2
Because we are interst by fit the homogeneity within this Experiment
factor.
We use this factor as random intercept.
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Experiment.RData")
# model.num.05 <- mixed(Mesure.TP ~ (Body + ns(Time,7))^2 +
# (1|Subject) + (ns(Time,7)|Body) + (1|Localization),
# family = Gamma("identity"),
# data = datas.raw.num,
# method = "LRT")
#
# save(model.num.05,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Experiment.RData")First we explore the difference between previous model.
anova(model.num.04, model.num.05)## Data: data
## Models:
## model.num.04: Mesure.TP ~ (Body + ns(Time, 7))^2 + (1 | Subject) + (ns(Time,
## model.num.04: 7) | Body)
## model.num.05: Mesure.TP ~ (Body + ns(Time, 7))^2 + (1 | Subject) + (ns(Time,
## model.num.05: 7) | Body) + (1 | Localization)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model.num.04 66 12683 13136 -6275.4 12551
## model.num.05 67 12685 13145 -6275.4 12551 0.0029 1 0.9572
There is no significative differences between these models. And we see an BIC and AIC increase with the last model. So we conclude to the lack of relevance for this new model.
Without surpise, we found the same significativity:
Anova(model.num.05$full_model, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Mesure.TP
## Chisq Df Pr(>Chisq)
## (Intercept) 8.4759 1 0.003599 **
## Body 1.5386 3 0.673392
## ns(Time, 7) 122.2503 6 < 2.2e-16 ***
## Body:ns(Time, 7) 54.0876 18 1.778e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.05$full_model, terms = c("Time", "Body")))em <- emmeans(model.num.05$full_model, ~ Time | Body,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")## Body = FOREARM:
## contrast estimate SE df z.ratio p.value
## 0_15 0.11904 0.0592 Inf 2.010 0.3551
## 0_25 0.16529 0.0542 Inf 3.048 0.0184
## 0_35 0.15415 0.0569 Inf 2.710 0.0539
## 0_45 0.10234 0.0599 Inf 1.709 0.6991
## 0_55 0.14000 0.0547 Inf 2.561 0.0835
## 0_65 0.16375 0.0571 Inf 2.867 0.0332
## 0_75 0.10408 0.0663 Inf 1.570 0.9305
## 0_85 0.18800 0.0627 Inf 3.000 0.0216
##
## Body = ARM:
## contrast estimate SE df z.ratio p.value
## 0_15 -0.05128 0.0721 Inf -0.711 1.0000
## 0_25 0.04127 0.0603 Inf 0.684 1.0000
## 0_35 0.12735 0.0621 Inf 2.051 0.3220
## 0_45 0.12308 0.0644 Inf 1.911 0.4478
## 0_55 0.16879 0.0572 Inf 2.949 0.0255
## 0_65 0.16176 0.0630 Inf 2.568 0.0819
## 0_75 0.03888 0.0722 Inf 0.538 1.0000
## 0_85 0.07571 0.0720 Inf 1.052 1.0000
##
## Body = CHEEK:
## contrast estimate SE df z.ratio p.value
## 0_15 -0.08016 0.0739 Inf -1.085 1.0000
## 0_25 0.03640 0.0635 Inf 0.573 1.0000
## 0_35 0.21999 0.0626 Inf 3.513 0.0035
## 0_45 0.35803 0.0588 Inf 6.085 <.0001
## 0_55 0.38905 0.0571 Inf 6.814 <.0001
## 0_65 0.29121 0.0648 Inf 4.493 0.0001
## 0_75 0.11289 0.0778 Inf 1.451 1.0000
## 0_85 0.05732 0.0826 Inf 0.694 1.0000
##
## Body = LIP:
## contrast estimate SE df z.ratio p.value
## 0_15 0.00204 0.0667 Inf 0.031 1.0000
## 0_25 0.00451 0.0595 Inf 0.076 1.0000
## 0_35 0.08774 0.0624 Inf 1.407 1.0000
## 0_45 0.27841 0.0567 Inf 4.909 <.0001
## 0_55 0.37806 0.0523 Inf 7.224 <.0001
## 0_65 0.32106 0.0560 Inf 5.734 <.0001
## 0_75 0.15246 0.0643 Inf 2.370 0.1424
## 0_85 0.07143 0.0702 Inf 1.017 1.0000
##
## P value adjustment: bonferroni method for 8 tests
Two separate analysis
Reviewers ask to analyse Experiment 1 and Experiment 2 datas
separatly. They are affraid that our GLMM hide an Experiment
effect.
So, we split (physically) datas according to Experiment 1 and Experiment
2 and we fit 2 models:
Experiment 1 (LIP & CHEEK)
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Supp.RData")
# model.num.06 <- mixed(Mesure.TP ~ (Body + ns(Time,7))^2 +
# (1|Subject) + (ns(Time,7)|Body),
# family = Gamma("identity"),
# data = filter(datas.raw.num, Localization == "Loc1"),
# method = "LRT")
Anova(model.num.06$full_model, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Mesure.TP
## Chisq Df Pr(>Chisq)
## (Intercept) 1.5181 1 0.2179
## Body 0.1043 1 0.7467
## ns(Time, 7) 69.3255 6 5.622e-13 ***
## Body:ns(Time, 7) 3.8732 6 0.6938
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.06$full_model, terms = c("Time", "Body")))em <- emmeans(model.num.06$full_model, "Time",
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")## contrast estimate SE df z.ratio p.value
## 0_15 -0.0333 0.0535 Inf -0.622 1.0000
## 0_25 0.0289 0.0466 Inf 0.620 1.0000
## 0_35 0.1628 0.0473 Inf 3.446 0.0046
## 0_45 0.3238 0.0441 Inf 7.349 <.0001
## 0_55 0.3811 0.0414 Inf 9.201 <.0001
## 0_65 0.3015 0.0461 Inf 6.536 <.0001
## 0_75 0.1377 0.0538 Inf 2.559 0.0839
## 0_85 0.0747 0.0580 Inf 1.287 1.0000
##
## Results are averaged over the levels of: Body
## P value adjustment: bonferroni method for 8 tests
The same results as before
Experiment 2 (FOREAM & ARM)
# model.num.07 <- mixed(Mesure.TP ~ (Body + ns(Time,7))^2 +
# (1|Subject) + (ns(Time,7)|Body),
# family = Gamma("identity"),
# data = filter(datas.raw.num, Localization == "Loc2"),
# method = "LRT")
Anova(model.num.07$full_model, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Mesure.TP
## Chisq Df Pr(>Chisq)
## (Intercept) 29.715 1 5.005e-08 ***
## Body 3054.248 1 < 2.2e-16 ***
## ns(Time, 7) 93350.704 6 < 2.2e-16 ***
## Body:ns(Time, 7) 23669.669 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.07$full_model, terms = c("Time", "Body")))em <- emmeans(model.num.07$full_model, ~ Time,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")## contrast estimate SE df z.ratio p.value
## 0_15 0.0291 0.0338 Inf 0.859 1.0000
## 0_25 0.1284 0.0215 Inf 5.962 <.0001
## 0_35 0.1383 0.0137 Inf 10.099 <.0001
## 0_45 0.1204 0.0131 Inf 9.194 <.0001
## 0_55 0.1602 0.0135 Inf 11.861 <.0001
## 0_65 0.1601 0.0135 Inf 11.881 <.0001
## 0_75 0.0749 0.0132 Inf 5.695 <.0001
## 0_85 0.1357 0.0111 Inf 12.262 <.0001
##
## Results are averaged over the levels of: Body
## P value adjustment: bonferroni method for 8 tests
# save(model.num.06, model.num.07,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Supp.RData")The same results as before exept for :
ARM 0_65 ms
FOREARM
0_55
becomes significant
Four separate analysis
The same as previously. Reviwers ask to analyse each body parts
independently.
So we split datas according to Body levels and we fit 4 models:
LIP
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Four_separate.RData")
# model.num.08 <- mixed(Mesure.TP ~ ns(Time,7) + (1|Subject),
# family = Gamma("identity"),
# data = filter(datas.raw.num, Body == "LIP"),
# method = "LRT")
Anova(model.num.08$full_model, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Mesure.TP
## Chisq Df Pr(>Chisq)
## (Intercept) 3.9684 1 0.04636 *
## ns(Time, 7) 69.1470 6 6.116e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.08$full_model, terms = c("Time")))em <- emmeans(model.num.08$full_model, ~ Time,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")## contrast estimate SE df z.ratio p.value
## 0_15 0.00165 0.0675 Inf 0.024 1.0000
## 0_25 0.00848 0.0598 Inf 0.142 1.0000
## 0_35 0.09688 0.0627 Inf 1.545 0.9795
## 0_45 0.28618 0.0579 Inf 4.939 <.0001
## 0_55 0.36836 0.0529 Inf 6.966 <.0001
## 0_65 0.30150 0.0578 Inf 5.213 <.0001
## 0_75 0.14755 0.0645 Inf 2.288 0.1772
## 0_85 0.07338 0.0709 Inf 1.035 1.0000
##
## P value adjustment: bonferroni method for 8 tests
The same results as before
CHEEK
# model.num.09 <- mixed(Mesure.TP ~ ns(Time,7) + (1|Subject),
# family = Gamma("identity"),
# data = filter(datas.raw.num, Body == "CHEEK"),
# method = "LRT")
Anova(model.num.09$full_model, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Mesure.TP
## Chisq Df Pr(>Chisq)
## (Intercept) 1.130 1 0.2878
## ns(Time, 7) 57.306 6 1.584e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.09$full_model, terms = c("Time")))em <- emmeans(model.num.09$full_model, ~ Time,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")## contrast estimate SE df z.ratio p.value
## 0_15 -0.0626 0.0853 Inf -0.734 1.0000
## 0_25 0.0544 0.0732 Inf 0.743 1.0000
## 0_35 0.2319 0.0720 Inf 3.222 0.0102
## 0_45 0.3628 0.0679 Inf 5.345 <.0001
## 0_55 0.3901 0.0658 Inf 5.927 <.0001
## 0_65 0.2943 0.0748 Inf 3.935 0.0007
## 0_75 0.1261 0.0895 Inf 1.409 1.0000
## 0_85 0.0901 0.0950 Inf 0.949 1.0000
##
## P value adjustment: bonferroni method for 8 tests
The same results as before
FOREARM
# model.num.10 <- mixed(Mesure.TP ~ ns(Time,7) + (1|Subject),
# family = Gamma("identity"),
# data = filter(datas.raw.num, Body == "FOREARM"),
# method = "LRT")
Anova(model.num.10$full_model, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Mesure.TP
## Chisq Df Pr(>Chisq)
## (Intercept) 0.4698 1 0.493084
## ns(Time, 7) 19.7598 6 0.003055 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.10$full_model, terms = c("Time")))em <- emmeans(model.num.10$full_model, ~ Time,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")## contrast estimate SE df z.ratio p.value
## 0_15 0.136 0.0525 Inf 2.597 0.0752
## 0_25 0.182 0.0485 Inf 3.758 0.0014
## 0_35 0.167 0.0509 Inf 3.276 0.0084
## 0_45 0.112 0.0535 Inf 2.094 0.2902
## 0_55 0.143 0.0485 Inf 2.945 0.0259
## 0_65 0.168 0.0507 Inf 3.305 0.0076
## 0_75 0.120 0.0586 Inf 2.053 0.3209
## 0_85 0.194 0.0555 Inf 3.503 0.0037
##
## P value adjustment: bonferroni method for 8 tests
The same results as before
Exept :
0_55 become
significant
ARM
# model.num.11 <- mixed(Mesure.TP ~ ns(Time,7) + (1|Subject),
# family = Gamma("identity"),
# data = filter(datas.raw.num, Body == "ARM"),
# method = "LRT")
Anova(model.num.11$full_model, type = "III")## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Mesure.TP
## Chisq Df Pr(>Chisq)
## (Intercept) 0.0054 1 0.94140
## ns(Time, 7) 16.7440 6 0.01027 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ggpredict(model.num.11$full_model, terms = c("Time")))em <- emmeans(model.num.11$full_model, ~ Time,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
contrast(em, des_cont, adjust = "bonferroni")## contrast estimate SE df z.ratio p.value
## 0_15 -0.0458 0.0732 Inf -0.626 1.0000
## 0_25 0.0530 0.0618 Inf 0.857 1.0000
## 0_35 0.1219 0.0546 Inf 2.232 0.2049
## 0_45 0.1582 0.0595 Inf 2.661 0.0623
## 0_55 0.1817 0.0571 Inf 3.185 0.0116
## 0_65 0.1383 0.0545 Inf 2.539 0.0889
## 0_75 0.0814 0.0606 Inf 1.345 1.0000
## 0_85 0.0717 0.0681 Inf 1.053 1.0000
##
## P value adjustment: bonferroni method for 8 tests
#
# save(model.num.08, model.num.09, model.num.10, model.num.11,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_Four_separate.RData")The same results as before
Power
Reviewers ask effect size or power estimation of our
statistics.
This is a tricky question.
We usetwo indirect ways to answer:
Rsquare
Rsquare give an idea of fit quality of our models. So we compare the Rsquare.
library(MuMIn)
body.name <- c("ARM", "FOREARM", "LIP", "CHEEK")
glmer.R2 <- c(r.squaredGLMM(glmer(Mesure.TP ~ ns(Time,7) + (1|Subject), data = filter(datas.raw.num, Body == "ARM")))[2],
r.squaredGLMM(glmer(Mesure.TP ~ ns(Time,7) + (1|Subject), data = filter(datas.raw.num, Body == "FOREARM")))[2],
r.squaredGLMM(glmer(Mesure.TP ~ ns(Time,7) + (1|Subject), data = filter(datas.raw.num, Body == "LIP")))[2],
r.squaredGLMM(glmer(Mesure.TP ~ ns(Time,7) + (1|Subject), data = filter(datas.raw.num, Body == "CHEEK")))[2])
aov.R2<-c(summary(lm(Mesure.TP ~ as.factor(Time), data = filter(datas.raw.num, Body == "ARM")))$r.squared,
summary(lm(Mesure.TP ~ as.factor(Time), data = filter(datas.raw.num, Body == "FOREARM")))$r.squared,
summary(lm(Mesure.TP ~ as.factor(Time), data = filter(datas.raw.num, Body == "LIP")))$r.squared,
summary(lm(Mesure.TP ~ as.factor(Time), data = filter(datas.raw.num, Body == "CHEEK")))$r.squared)
data.frame(body = body.name,
glmer = glmer.R2,
aov = aov.R2) %>%
mutate(gain = (glmer-aov)/aov)## body glmer aov gain
## 1 ARM 0.04769246 0.010683196 3.464250
## 2 FOREARM 0.04216606 0.006509172 5.477945
## 3 LIP 0.05642969 0.020136078 1.802417
## 4 CHEEK 0.05296593 0.016857169 2.142041
data.frame(body = c("LIP & CHEEK", "ARM & FOREARM"),
glmer = c(r.squaredGLMM(glmer(Mesure.TP ~ (Body + ns(Time,7))^2 + (1|Subject) + (ns(Time,7)|Body),data = filter(datas.raw.num, Localization == "Loc1")))[2],
r.squaredGLMM(glmer(Mesure.TP ~ (Body + ns(Time,7))^2 + (1|Subject) + (ns(Time,7)|Body),data = filter(datas.raw.num, Localization == "Loc2")))[2]),
aov = c(summary(lm(Mesure.TP ~ (Body + as.factor(Time))^2 ,data = filter(datas.raw.num, Localization == "Loc1")))$r.squared,
summary(lm(Mesure.TP ~ (Body + as.factor(Time))^2 ,data = filter(datas.raw.num, Localization == "Loc2"))))$r.squared) %>%
mutate(gain = (glmer-aov)/aov)## body glmer aov gain
## 1 LIP & CHEEK 0.6152456 0.008596999 70.56516
## 2 ARM & FOREARM 0.6094596 0.008596999 69.89213
data.frame(body = c("LIP & CHEEK & ARM & FOREARM"),
glmer = r.squaredGLMM(glmer(Mesure.TP ~ (Body + ns(Time,7))^2 + (1|Subject) + (ns(Time,7)|Body) + (1|Localization), data = datas.raw.num))[2],
aov = summary(lm(Mesure.TP ~ (Body + as.factor(Time))^2, data = datas.raw.num))$r.squared) %>%
mutate(gain = (glmer-aov)/aov)## body glmer aov gain
## 1 LIP & CHEEK & ARM & FOREARM 0.7248271 0.01420603 50.02248
The more we split datas, the less models fit with datas.
Bayes framwork
Another way to check the models quality is to use the same models with bayesian estimation of fit. These model give us some metrics to compare the models quality.
1 model, 4 body parts
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_4.RData")
# model.stan.00 <- stan_glmer(Mesure.TP ~ (Body + ns(Time,3))^2 + (1|Subject) + (ns(Time,7)|Body) + (1|Localization), family = Gamma("identity"), data = datas.raw.num,
# diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
#
# save(model.stan.00,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_4.RData")
plot(ggpredict(model.stan.00, terms = c("Time", "Body"))) +
ylim(c(0.6, 1.25))emm <- emmeans(model.stan.00, ~ Time | Body,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 80)))
BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.00)
as.data.frame(p_direction(contrast(emm, des_cont))) %>%
mutate(BF = exp(BF.emm$log_BF)) %>%
mutate(BayesSign = case_when(BF > 100 ~ "***",
BF > 10 & BF < 100 ~ "**",
BF > 3 & BF < 10 ~ "*",
BF > 1 & BF < 3 ~ ".",
BF < 1 ~ "")) -> Tests.Pairs.Results
datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )rm(Tests.Pairs.Results, emm, BF.emm)2 models, 2 body parts each
LIP & CHEEK
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_2_LIP_CHEEK.RData")
# model.stan.01 <- stan_glmer(Mesure.TP ~ (Body + ns(Time,3))^2 + (1|Subject) + (ns(Time,7)|Body), family = Gamma("identity"),
# data = filter(datas.raw.num, Localization == "Loc1"),
# diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
#
# save(model.stan.01,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_2_LIP_CHEEK.RData")
plot(ggpredict(model.stan.01, terms = c("Time", "Body")))emm <- emmeans(model.stan.01, ~ Time | Body,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.01)
as.data.frame(p_direction(contrast(emm, des_cont))) %>%
mutate(BF = exp(BF.emm$log_BF)) %>%
mutate(BayesSign = case_when(BF > 100 ~ "***",
BF > 10 & BF < 100 ~ "**",
BF > 3 & BF < 10 ~ "*",
BF > 1 & BF < 3 ~ ".",
BF < 1 ~ "")) -> Tests.Pairs.Results
datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )rm(Tests.Pairs.Results, emm, BF.emm)FOREARM & ARM
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_2_FOREARM_ARM.RData")
# model.stan.02 <- stan_glmer(Mesure.TP ~ (Body + ns(Time,3))^2 + (1|Subject) + (ns(Time,7)|Body), family = Gamma("identity"),
# data = filter(datas.raw.num, Localization == "Loc2"),
# diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
# save(model.stan.02,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_2_FOREARM_ARM.RData")
plot(ggpredict(model.stan.02, terms = c("Time", "Body")))emm <- emmeans(model.stan.02, ~ Time | Body,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.02)
as.data.frame(p_direction(contrast(emm, des_cont))) %>%
mutate(BF = exp(BF.emm$log_BF)) %>%
mutate(BayesSign = case_when(BF > 100 ~ "***",
BF > 10 & BF < 100 ~ "**",
BF > 3 & BF < 10 ~ "*",
BF > 1 & BF < 3 ~ ".",
BF < 1 ~ "")) -> Tests.Pairs.Results
datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )rm(Tests.Pairs.Results, emm, BF.emm)4 models
LIP
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_LIP.RData")
# model.stan.03 <- stan_glmer(Mesure.TP ~ ns(Time,3) + (1|Subject), family = Gamma("identity"),
# data = filter(datas.raw.num, Body == "LIP"),
# diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
#
#
#
# save(model.stan.03,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_LIP.RData")
plot(ggpredict(model.stan.03, terms = c("Time")))emm <- emmeans(model.stan.03, ~ Time,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.03)
as.data.frame(p_direction(contrast(emm, des_cont))) %>%
mutate(BF = exp(BF.emm$log_BF)) %>%
mutate(BayesSign = case_when(BF > 100 ~ "***",
BF > 10 & BF < 100 ~ "**",
BF > 3 & BF < 10 ~ "*",
BF > 1 & BF < 3 ~ ".",
BF < 1 ~ "")) -> Tests.Pairs.Results
datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )rm(Tests.Pairs.Results, emm, BF.emm)CHEEK
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_CHEEK.RData")
# model.stan.04 <- stan_glmer(Mesure.TP ~ ns(Time,3) + (1|Subject), family = Gamma("identity"),
# data = filter(datas.raw.num, Body == "CHEEK"),
# diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
#
# save(model.stan.04,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_CHEEK.RData")
plot(ggpredict(model.stan.04, terms = c("Time")))emm <- emmeans(model.stan.04, ~ Time,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.04)
as.data.frame(p_direction(contrast(emm, des_cont))) %>%
mutate(BF = exp(BF.emm$log_BF)) %>%
mutate(BayesSign = case_when(BF > 100 ~ "***",
BF > 10 & BF < 100 ~ "**",
BF > 3 & BF < 10 ~ "*",
BF > 1 & BF < 3 ~ ".",
BF < 1 ~ "")) -> Tests.Pairs.Results
datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )rm(Tests.Pairs.Results, emm, BF.emm)FOREARM
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_FOREARM.RData")
# model.stan.05 <- stan_glmer(Mesure.TP ~ ns(Time,3) + (1|Subject), family = Gamma("identity"),
# data = filter(datas.raw.num, Body == "FOREARM"),
# diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
# save(model.stan.05,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_FOREARM.RData")
plot(ggpredict(model.stan.05, terms = c("Time")))emm <- emmeans(model.stan.05, ~ Time,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.05)
as.data.frame(p_direction(contrast(emm, des_cont))) %>%
mutate(BF = exp(BF.emm$log_BF)) %>%
mutate(BayesSign = case_when(BF > 100 ~ "***",
BF > 10 & BF < 100 ~ "**",
BF > 3 & BF < 10 ~ "*",
BF > 1 & BF < 3 ~ ".",
BF < 1 ~ "")) -> Tests.Pairs.Results
datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )rm(Tests.Pairs.Results, emm, BF.emm)ARM
load(file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_ARM.RData")
# model.stan.06 <- stan_glmer(Mesure.TP ~ ns(Time,3) + (1|Subject), family = Gamma("identity"),
# data = filter(datas.raw.num, Body == "ARM"),
# diagnostic_file = "/Users/romain/Desktop/diagnostic_00.csv")
#
# save(model.stan.06,
# file = "/Users/romain/Study/Reilly_Karen/Face_hand_sensori/Scripts/R/save_models/Models_Time_Nume_stan_1_ARM.RData")
plot(ggpredict(model.stan.06, terms = c("Time")))emm <- emmeans(model.stan.06, ~ Time,
at = list(Time = c(0, 15, 25, 35, 45, 55, 65, 75, 85)))
BF.emm <- bayesfactor_parameters(contrast(emm, des_cont), prior = model.stan.06)
as.data.frame(p_direction(contrast(emm, des_cont))) %>%
mutate(BF = exp(BF.emm$log_BF)) %>%
mutate(BayesSign = case_when(BF > 100 ~ "***",
BF > 10 & BF < 100 ~ "**",
BF > 3 & BF < 10 ~ "*",
BF > 1 & BF < 3 ~ ".",
BF < 1 ~ "")) -> Tests.Pairs.Results
datatable(Tests.Pairs.Results, filter="top", options = list(pageLength = 12, scrollX=T) )rm(Tests.Pairs.Results, emm, BF.emm)Version
Here, all packages en R version used on this analysis.
print(sessionInfo())## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
##
## attached base packages:
## [1] splines stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] MuMIn_1.43.17 bayestestR_0.12.1
## [3] bayesplot_1.9.0 brms_2.17.0
## [5] rstanarm_2.21.1 Rcpp_1.0.6
## [7] sjmisc_2.8.9 sjPlot_2.8.10
## [9] glmmTMB_1.1.3 DHARMa_0.4.4
## [11] ggeffects_1.1.1 ggpubr_0.4.0
## [13] emmeans_1.7.3 lmerTest_3.1-3
## [15] multcomp_1.4-16 TH.data_1.0-10
## [17] MASS_7.3-53 survival_3.2-7
## [19] mvtnorm_1.1-1 phia_0.2-1
## [21] car_3.0-10 carData_3.0-4
## [23] LMERConvenienceFunctions_3.0 RVAideMemoire_0.9-81-2
## [25] afex_1.1-1 lme4_1.1-26
## [27] Matrix_1.2-18 forcats_0.5.1
## [29] stringr_1.4.0 dplyr_1.0.8
## [31] purrr_0.3.4 readr_1.4.0
## [33] tidyr_1.2.0 tibble_3.1.6
## [35] tidyverse_1.3.1 ggplot2_3.3.6
## [37] openxlsx_4.2.3 readxl_1.3.1
## [39] readtext_0.80 DT_0.17
## [41] epuRate_0.1 rmarkdown_2.11
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.1 tidyselect_1.1.2 htmlwidgets_1.5.3
## [4] grid_4.0.3 munsell_0.5.0 codetools_0.2-16
## [7] effectsize_0.4.4 statmod_1.4.35 miniUI_0.1.1.1
## [10] withr_2.4.1 Brobdingnag_1.2-6 colorspace_2.0-0
## [13] highr_0.8 knitr_1.39 rstudioapi_0.13
## [16] stats4_4.0.3 ggsignif_0.6.1 labeling_0.4.2
## [19] LCFdata_2.0 rstan_2.21.1 farver_2.1.0
## [22] datawizard_0.4.0 bridgesampling_1.1-2 coda_0.19-4
## [25] vctrs_0.4.0 generics_0.1.0 xfun_0.30
## [28] R6_2.5.0 markdown_1.1 rmdformats_1.0.2
## [31] fields_11.6 gamm4_0.2-6 projpred_2.0.2
## [34] logspline_2.1.16 assertthat_0.2.1 promises_1.2.0.1
## [37] scales_1.1.1 gtable_0.3.0 processx_3.5.0
## [40] spam_2.6-0 sandwich_3.0-0 rlang_1.0.2
## [43] Rmisc_1.5 rstatix_0.7.0 TMB_1.7.19
## [46] checkmate_2.0.0 broom_0.8.0 inline_0.3.19
## [49] yaml_2.2.1 reshape2_1.4.4 abind_1.4-5
## [52] modelr_0.1.8 gtsummary_1.6.0 threejs_0.3.3
## [55] crosstalk_1.1.1 backports_1.2.1 httpuv_1.5.5
## [58] rsconnect_0.8.24 tensorA_0.36.2 tools_4.0.3
## [61] bookdown_0.22 ellipsis_0.3.2 RColorBrewer_1.1-2
## [64] posterior_1.1.0 jquerylib_0.1.3 ggridges_0.5.3
## [67] plyr_1.8.6 base64enc_0.1-3 ps_1.6.0
## [70] prettyunits_1.1.1 zoo_1.8-9 haven_2.4.3
## [73] fs_1.5.2 magrittr_2.0.3 data.table_1.14.0
## [76] colourpicker_1.1.1 reprex_2.0.1 matrixStats_0.58.0
## [79] hms_1.0.0 shinyjs_2.0.0 mime_0.10
## [82] evaluate_0.15 xtable_1.8-4 shinystan_2.5.0
## [85] rio_0.5.26 sjstats_0.18.1 gridExtra_2.3
## [88] rstantools_2.2.0 compiler_4.0.3 maps_3.3.0
## [91] gt_0.5.0 V8_3.4.2 crayon_1.4.1
## [94] minqa_1.2.4 StanHeaders_2.21.0-7 htmltools_0.5.2
## [97] mgcv_1.8-33 later_1.1.0.1 RcppParallel_5.1.4
## [100] lubridate_1.7.10 DBI_1.1.1 sjlabelled_1.1.7
## [103] dbplyr_2.1.1 broom.helpers_1.7.0 boot_1.3-25
## [106] cli_3.2.0 parallel_4.0.3 insight_0.17.0
## [109] dotCall64_1.0-1 igraph_1.2.6 pkgconfig_2.0.3
## [112] numDeriv_2016.8-1.1 foreign_0.8-80 xml2_1.3.2
## [115] dygraphs_1.1.1.6 bslib_0.2.4 estimability_1.3
## [118] rvest_1.0.0 distributional_0.2.2 callr_3.7.0
## [121] digest_0.6.27 parameters_0.14.0 cellranger_1.1.0
## [124] curl_4.3 commonmark_1.7 shiny_1.6.0
## [127] gtools_3.8.2 nloptr_1.2.2.2 lifecycle_1.0.1
## [130] nlme_3.1-149 jsonlite_1.7.2 fansi_0.4.2
## [133] pillar_1.7.0 lattice_0.20-41 loo_2.4.1
## [136] fastmap_1.1.0 httr_1.4.2 pkgbuild_1.2.0
## [139] glue_1.6.2 xts_0.12.1 zip_2.1.1
## [142] shinythemes_1.2.0 stringi_1.5.3 sass_0.4.1
## [145] performance_0.9.0